Numerical Computation of First-Passage Times of Increasing Levy 

Processes 



Mark Veillette and Murad S. Taqqu 
Boston University 

April 27, 2009 



Abstract 

Let {D(s), s > 0} be a non-decreasing Levy process. The first-hitting time process {E(t) t > 0} 
(which is sometimes referred to as an inverse subordinator) defined by E(t) = inf{s : D(s) > t} is 
a process which has arisen in many applications. Of particular interest is the mean first-hitting time 
U(t) = E.E(t). This function characterizes all finite-dimensional distributions of the process E. The 
function U can be calculated by inverting the Laplace transform of the function (7(A) = (A0(A)) _1 , 
where <j> is the Levy exponent of the subordinator D. In this paper, we give two methods for computing 
numerically the inverse of this Laplace transform. The first is based on the Bromwich integral and the 
second is based on the Post-Widder inversion formula. The software written to support this work is 
available from the authors and we illustrate its use at the end of the paper. 



1 Introduction 

Let {D{s), s > 0} be a Levy subordinator, that is, a non-decreasing Levy process starting from 0, which is 
continuous from the right with left limits. This process has stationary and independent increments and is 
characterized by its Laplace Transform 

Ee -AD(«) = e -«*(A) j A > _ ^) 

The function <f> above is known as the Levy exponent (or Laplace exponent) and is given by the Levy- 
Khintchinc formula: 

<KA)=/iA+/ (l-e- Xx )U(dx), (2) 

J(0,oo) 

where fj, > is the drift and the Levy measure II is a measure on R + U{0} which satisfies / °°(1 /\x)H(dx) < oo 
(see [3], [9] or [10]). 
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Consider the inverse subordinator {E(t), t > 0}, which is given by the first-passage time of D: 

E(t) = inf{s : D(s) > t}. (3) 

The process E(t) is non-decreasing, and its sample paths are a.s. continuous if and only if D is strictly 
increasing. Also, E is, in general, non-Markovian with non-stationary and non-independent increments. 

Inverse subordinators are of particular interest in the study of fractional kinetics and the scaling limits 
of continuous time random walks, [7J, [B], [2D], [23], [?]. Here, a Markov process, {X(t),t > 0}, is time- 
changed with an inverse subordinator, yielding the new process M(t) = X(E(t)), t > 0. For certain choices 
of subordinators D, using E as a time change gives a model of anomalous diffusion, where the mean-squared 
displacement of M grows non-linearly in time. For example, if D is an a-stable subordinator, that is, the 
subordinator whose Levy exponent is given by 

0(A) = A", < a < 1, (4) 

then the mean of E(i) is given by the power law KE(t) ~ t a . Other types of non-linear behavior are 
also possible. In [20j . Meerschaert considers the subordinator with Levy exponent given by the following 
generalization of (g|): 

0(A) - / X dp((3), (5) 

where p is a probability measure on (0, 1). Depending on the choice of p, the mean of E in this case can 
grow at logarithmic rates. Meerschaert also shows in this case that the transition density of the process M 
solves a time-fractional diffusion equation whose order is distributed according to p. 

Inverse Subordinators have also appeared in many other areas of probability theory. Early work regarding 
the joint distribution of E(t) and D(E{t)) was done in [15] and [19] . More recently, in [18], Kaj and Martin- 
Lof show that a scaled sum of these processes converges weakly to another non-stable and non-Gaussian 
process. An application of inverse subordinators in modeling of foreign exchange markets was considered in 
[27] . Distributional properties of inverse subordinators were studied in [21j , which drew upon a connection 
with Cox processes. A different approach using differential equations to characterize the joint distribution 
function of the process E was used in [2D] and [7J . 

An important function in the study of inverse subordinators is the so-called renewal function, U(t), t > 0, 
which is given by the mean of the inverse subordinator, U(t) = ME(t). It has be shown by different methods 
that this function characterizes all finite-dimensional distributions of the process E. The Laplace transform 
of U is given simply in terms of the Levy exponent 0, however, inverting this Laplace transform to obtain 
U in closed form is not always possible. In this paper, we provide two numerical methods for obtaining 
this mean first-passage time which take as input the drift [i and the Levy measure II of the corresponding 
subordinator D. Once U is obtained, higher order moments can be calculated by methods described in [26j 
or [21]. 

This paper is organized as follows: We begin in Section [2] by giving a brief background on the theory 
of inverse subordinators, as well as describe the importance of the function U. In Section [3] we develop 
two methods for calculating U numerically and test them in cases where U can be computed analytically. 
In Section 0] we apply these techniques to the following examples: Poisson process with drift (Section l4.1[) . 
Compound Poisson process with Pareto jumps ( Section 14. 2p . special cases of the "mixed" a-stable process 
( Sect ion |4~3]) . and the generalized inverse Gaussian Levy process fSection l4~4]) . For each example, we compute 
one and two-time moments of E. Methods to compute U for each example are implemented in MATLAB. 
The software for doing this is freely available from the authors and its use is described in Section [6J 
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2 Background 

We start by recalling the fundamental relationship between a subordinates and its inverse. If D is strictly 
increasing, then for t\, t%, • . • , t n and Si, S2, ■ ■ ■ , s n positive, 

{D( Sl ) <U,i = l,...n} = {E(U) >s h i=l,..., n}. (6) 

If D is not strictly increasing, then the above relationship holds off a set of measure in (ti, t%, . . . , t n ) (see 

1251). 

As in the introduction, define the renewal function to be the mean of the inverse subordinator, U (t) = 
EE(t) fori > 0. Also, let H s (t) = P[D(s) < t]. From the Levy-Khintchine formula, we have J °° e~ xt dH s (t) = 
e -s<t>W _ This together with ([6]) lets us compute the Laplace transform, U, of U: 



U(\) = I U(t) e - Xt dt (7) 

/ P[E(t) > s]e- xt dsdt (8) 
Jo 

/ P[D{s) < t]er xt dsdt (9) 
Jo 

/ -e- xt dH s (t)ds (10) 

JO * 



Thus, U characterizes the process E (since <j> characterizes D). Since E is non-decreasing, we can define 
the Borel measure dU which is induced by U, which is commonly referred to as the renewal measure. Notice 
that the renewal measure has the following property. For a.e. < a < b, we have 

l M (r)dU(r) = U(b)-U(a) 

= / P[E(b) > s) - P[E(a) > s]ds 
Jo 

poo 

= / P[D{s) <b}- P[D(s) < a]ds 
Jo 

poo 

= E (l ( _ 00!6] ( J D(s))-l ( _ 00!0] (£>(s)))ds 
Jo 

poo 

= E / l {aM (D(s))ds. (12) 



By approximating with step functions, this relationship can be extended to continuous functions g as 
Jo°° 9( T )dU(T) — E Jp 00 g(D(s))ds. With this, we see the Laplace transform of dU is given by 



/•oo />oo />OC 

/ e~ TX dU(T) =E e~ XD ^ds = / e'^ds 
Jo Jo Jo 



0(A)' 



(13) 



The pair U and dU can be used to compute all joint moments of the process E. For non-negative integers 
mi , . . . , m„ define 

U(t u . . . ,t„;m 1 , . . . ,m„) = EE(h) mi E(t 2 ) m2 . . . E(t n ) m " . (14) 
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The order of the moment in (fTT)) is denned to be N = Y^i=i m i- The following theorem from [55] gives the 
n-time Laplace transform U{\\, . . . , A„; mi, . . . , m n ) in terms of a strictly lower older moment. 

Theorem 2.1 Let D be a general Levy subordinator with Levy exponent <f> and let E be the inverse subor- 
dinator of D. The n-time Laplace Transform of the N th order moment U(ti, . . . ,t n ; mi, . . . , m„) defined in 
{1$ is given by 



U(Xi, . .., X n ;mi, ... , m„) = — — ■ — — m,iU(\i, ...,X n ;mi,... ,m l ^ 1 ,m l - l,m i+1 , . . . , m„). 

0(Ai H h A„) ^ 

(15) 

Notice that the Laplace transform of U{t\,. .. ,t n ;mi, . . . ,m„) is given as the product of 1/0 and the 
Laplace transform of a strictly lower order moment. Taking inverse Laplace transforms, the N th order 
moment is given as the sum of convolutions 

71 />miiii ti 

U(ti,...,t n ;mi,... ,m„) = y^mj / f/(ti - r, . . . ,t n - r;rni, . . .,m i ^ 1 ,m l ~ l,m i+ i, . . .,m n )dU{r). 

<=i Jo 

(16) 

For full details, see 

Thus, if one knows the function U(t; 1) = U(t), t > 0, then all higher order moments can be obtained 
inductively. For example, the covariance is given by 

Cov(E(t 1 ),E(t 2 ))= {V(ti-T) + Ufa-T))dU(T)-U(ti)Ufa). (17) 



JO 

It is not always easy, however, to invert the Laplace transform (fTTj) to obtain U(t) analytically. In 
the following, we give two methods for numerically inverting this Laplace Transform given only the drift 
H and Levy measure LT of the subordinator D. Once U is obtained, the density of the renewal measure 
can be obtained by numerical differentiation and integral expressions such as (|17|) can be approximated by 
numerical integration. The first method is based on the Bromwhich integral which expresses the inverse 
Laplace transform as a path integral in the complex plane, and the second is based on the Post-Widdcr 
inversion formula, which expresses U as a limit of terms involving derivatives of U. 



3 Computing U(t) 



In Section [5] we saw that all moments of an inverse subordinator can be computed if one has first computed 
the renewal function U(t) = KE(t), which is given by the inverse Laplace transform of (fTTj) . In some cases, 
an analytical expression for U can be found, and in most cases, the asymptotics of U can be studied using a 
Tauberian Theorem. In this section we give two methods for calculating U numerically. The first is simple 
and precise, but is difficult to use when is a complicated function. The second is more robust, but requires 
smoothness in U to be effective. 

3.1 Method 1: Numerical Integration 

The first method involves calculating the inverse Laplace transform of U by numerically approximating the 
Bromwich integral: 



1 



2iri 



U{t) = — / e zt U{z)dz, (18) 



b-i 
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where b is chosen such that U is analytic in the region Rc(z) > b. Using the fact that U(t) = for t < 0, 
and symmetry properties of analytic functions which arc real valued on the real axis, (|18p simplifies into two 
equivalent expressions: 

2gbt poo —2e bt r °° 



U(t) = —j (Re(U(b + iu)) cos(ui)) du = — — J (lm(U(b + iu)) sin(ui)) du. (19) 

For details see [T]. 

From the Levy Khintchine formula @, we have 

cf>(b + iu) = fib + ifiu+ f ( 1 - e - x ( b+m A U(dx) (20) 



= (l J ' b + J q {1 - e~ xb cos(xu)) n(dx)^j +i \^iu- e~ xb sin(xu)Il(dx)\ (21) 
ee <f> r (b + iu) + i<j>i(b + iu). (22) 
Now, since U(X) = (A0(A)) _1 , a simple calculation gives 

(b z + u z )((p r (b + iu) z + <pi{b + iu) z ) 

~ . b4> t {b + iu) + u(f> r (b + iu) 

Im(U(b + i u)) = {b2 + u2){Mb + lu)2 + Mb + lury (24) 

Since 0(0) = and <f> is increasing in A, U(X) = (A^(A)) -1 has a singularity the origin. Therefore, we 
much choose b > above. Then, given /x and II, we evaluate the integrals in (|2ip to obtain <j) r and 4>i and 
use (|23|) and p4|) to obtain Re(U(b + iu)) and lm(U(b + iu)) for fixed u and then compute cither integral 
in (fl9|) to get U(t). Alternatively, if 4> is known in closed form, then (f> r and 4>i can be computed directly. 
This method works fairly well when 4> r and <j>i are easy to compute, for instance, in the case of the Poisson 
process with drift. The main problem with the method is when the integrands in (flT)]) and/or (|2~lj) are highly 
oscillatory and slowly decaying, causing most integration algorithms to converge very slowly. 

3.2 Method 2: Post-Widder Inversion 

The following method is based on the Post-Widder inversion formula (pQ, Theorem 2 or [T^], section VII. 6). 
In order to justify using this formula in our case, we state and prove this result under slightly weaker 
conditions than those found in [12] , 

Theorem 3.1 (Post-Widder Inversion) . Let u : R + — > K. + be a continuous function such that u(x)/x < C 
as x — > oo for some C > 0, and let u(X) be the Laplace transform of u. Then for t > 0, we have 

«"-JS^0)'^(f). 

where u^ denotes the k th derivative, k G Z + . 
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Proof. Fix t > 0. Let JQ, i = 1,2,... be iid gamma random variables with mean t and variance t 
and let M k — j Ei=i -^r The idea is to approximate w(t) by Ew(Affc). Notice that has a gamma 
distribution with mean t and variance t 2 /k. By the law of large numbers, M k — > d t as k — ► oo and by 
the continuous mapping theorem ([15], Theorem 5.10.4), u(Affc) u(t). Observe that we also have the 
following convergence in mean: 

Eu{M k ) -> u(t), (26) 

which follows from the uniform integrability of u(Mfc). This can be seen using the density of a gamma 
random variable and the assumption that u(x) < Cx as x — > oo. Indeed, we have 

/>oo 

lim supEu(M fc )l Mfc >a = hm sup / u(x)x k ~ 1 e~ kx/t dx (27) 

Q ->°° fe>l a-* 00 fe>l 7 a 

/•oo 

< lim supC / x k e~ kx/t dx (28) 



a — >oo 



fe>l 



lim supC / (xe~ x/t ) k dx (29) 
a ^°° fc>i 



lim C / (k" 1 /')^ (30) 
0. (31) 



Now, (|26|) implies 

i / h \ ^ r°° 

u{t) = hm Ew(Mfc) = lim — — — ( -J y u^z^V^da: (32) 



(-l)^ 1 (J£ 

k^L {k-l)l \t j 



lim ( 1 )' 8 ' C ^ " u^ik/t). (34) 



This verifies 



Remarks. 1. In order to justify using the Post-Widder formula, we must be sure that the renewal function 
U satisfies U(t)/t < C as t — > oo. Fortunately, the Renewal Theorem (see for instance, Proposition 3 in [2Tj ) 
implies that if D has finite mean, then U{t) will grow as 0{t) as t — > oo. If D has infinite mean, its mean 
first hitting time EE(t) will, on average, grow slower than in the finite mean case for large t. To see this, any 
subordinator D with infinite mean can be bounded below by a subordinator D with finite mean by truncating 
the large jumps of L0. The inverse, E, of D will be a.s. greater than E, and thus U(t) < U(t) = 0(t). 

2. One must check that U is continuous before using the Post-Widder formula. Proposition A.l in [26] 
implies that U is continuous if D is strictly increasing. 

Thus, calculating U(t) involves two steps: 
1. Calculating 



lr That is, D(s) < D(s) for all s > a.s. 
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for some set of integers ki, fe, • • • , k n (the choice of the fej's will be addressed later). 

2. Using the values of Uk t (t), i = 1, .. . n, to approximate the limit in (j2"S"l) . 

We will first focus on step 1 above. The main difficulty with using this method resides in the fact that 
(j2"5)) involves derivatives of arbitrarily high orders, which are often difficult to compute. Fortunately, in our 
case, this calculation can be done in a reasonable fashion, which we now outline. Recall Leibnitz's formula, 
which states that for smooth functions /, g, 



^/(*)fl(*) = E(J)f (M) (*)/ (0 (*)- 



(36) 



Let k > l,i > 0, and ^>(A) = l/4>(\). Applying Leibnitz's formula to U(X) = A 1- 0(A) with A = k/t gives 



(fc-l)! \t 



(k-l)\ \t 



k fc-1 



i=0 



fc-l 



i=0 



where the vectors Vfc, € R fc are given by 

(v fc ) t = 
(Wjfe)i = 



(-1)^ 



To compute the components of Wfc, we use an idea from 
function f = 0(A)^(A), we see that for any A > 0, 



E(V)( '-'^ifc 1 -' ) V <»/ 



i = 0, 1 
i = 0,...,&- 1. 



(38) 
(39) 

(40) 
(41) 



Using Leibnitz's formula on the unit 



E( •)^" i) (A)^ i) (A) = ^J0(A)V>(A) 



0» 



9A J 
dXi 

fi i = o, 
lo i>i 



i = o, 



i. 



(42) 

(43) 
(44) 



From this, we obtain the following matrix equation: 



/ 0(A) 
0'(A) 
0"(A) 

r(A) 





0(A) 
2^' (A) 
30" (A) 







0(A) 

30'(A) 0(A) 



V ^-^(A) (fc- l)0( fe - 2 )(A) 



\ 



0"(A) 
ip"'(X) 



( 1 \ 








0(A) ) V V' (fe - 1} (A) / V ) 



(45) 
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Choosing A = k/t, we conclude that satisfies the matrix equation GfcWfe 
K fc and 



0. 



1 < j < i < k 



e k , where e k = (1,0, ... ,0)' £ 

(46) 



Observe that the entries in G k are easily expressed in terms of the drift fj, and the Levy measure II corre- 
sponding to the subordinate!' D. Indeed, we have 



<j)^{k/t) = 



k 



1-e 



-kx jt 



U(dx), j = 0, 



-kx/t 



U(dx), j = l, 



(47) 



(-1) 



3 + 1 



i' —kx/t 

x J e ' 



IL(dx), j>2 



Thus, to compute wj, one needs to compute the entries of G k using (|46[) and (|47p . and then solve the matrix 
equation GfeW^ = e k . Notice that the integrals in (|47|) are much easier to compute than those in Section [3T 
since the integrands do not oscillate; they decay exponentially and are positive. Observe that Vfe, and hence 
Vfe • Wfc are easy to compute. Using this method, we get Uk^t) = Vfe ; ■ Wfe i in |35|) . 



Remark. The density U'(t) of the renewal measure is also of interest, for example in (|T7|) . It can be 
approximated as U is in (|35|) . since the Laplace transform of U'(t) is ^(A). One needs to compute f/;^ -1 ^ (k/t), 
which is obtained with no extra cost from (|45|) since it is the last component of w^. 



For step 2 above, we refer to the technique explained in [13] . To summarize, it has been shown ([2].|17j) 
that if U is smooth, then we have the following series expansion: 



u k (t) = u(t)+J2 



a m (t) 



(48) 



m— 1 



where a m (i) are remainder terms. Write hi = l/ki and let Unfa) = Uk^)- With this, 

oo 



implies 



(49) 



m— 1 



Our goal is to compute Uo(t) = U(t) = lim^oo Uk(t) with t fixed. To do so, we consider Uh(t) as a function 
of h and, given hi, hi, . . . , h n and Uh 1 , ■ ■ ■ , Uj ln , we write down the so-called "Lagrange polynomial" P n (h): 
this is the polynomial of degree n — 1 which passes through the n points (hi, Uhi (t)), i= 1,2, ... ,n, 



Pn(h)=J2[U h ~ hj 



i=l \jjti 



hi 



u hi (t)- 



(50) 



Observe that if h = h k for some k = 1,2, ... ,n, then all the summands in (|50|) vanish except for the term 
where i = k, which is equal to U) lk (t). Thus the polynomial P n (h) passes through the n points (hi, Uhi(t)), 
i = 1,2, ... ,n. 
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Observe that we cannot compute Uo(t) since h = corresponds to k — oo. However, with the method 
described above, we can compute Uh t for hi ^ and then approximate Uo(t) by P„(0). Since Uo(t) = U(t), 
U is then approximated by the linear combination 

n 

U(t)^P n (0)=J2c^U k M (51) 

where the weights c$ , calculated by setting h = in (|50|) . depend on the choice of {/ii}i=i@- 

Because the above method allows us to calculate U( k > for large k, we shall take ki = 2 I_1 , i = 1, 2, . . . , n. 
The corresponding weights are given by 



c 



(») 



n;=i(2-'"-i)n;=i 4 (2-'"-i)^ 



i = l,2, ...,n. (52) 



Thus, given a desired level of accuracy e > and t > fixed, we compute the sequence Pi(0),P2(0), . . . 
until |P„(0) — Pi-i(0)| < e and then set P(i) = P n (0). Typically this method converges with n < 9 for 
reasonable values of e. Since fc 10 = 512, it is often impossible to compute for i > 10 in double precision 
arithmetic, as noted in the remarks below. 

Note that for (|4"5|) to hold and for this method to be most effective, U must be sufficiently smooth. We 
found that for points at which U is not diffcrentiablc, the sequence P; had a slower rate of convergence. 

We close this section with two remarks regarding the computation of Uk{t). 

Remark. Caution should be taken when computing terms like k 3 / j\ for k,j large, since numbers like 
100 100 and 200! arc outside the range of double precision arithmetic. Instead, one should use expressions 
like k' ] I j\ = cxp(j log(fc) — J2l=i l°g(*))> because while k 3 and j\ may be large, their ratio may be small. 
Nevertheless, for large enough k and j (relation (|40p . for example, which requires computing k k ~ 1 /(k— 1)!), 
even the ratio may be too large. This happens for instance if k = feu = 2 10 = 1024, in which case 
1024 1023 /1023! « 6.5 x 10 442 . To be safe, we used at most k 10 = 512. 

Remark. A similar overflow problem can occur for k/t large. To avoid this, one can make the following 
adjustment to (f2"5|) . Let c > 0, then 

<-i>'-'mV-»m . (53) 



(k-l)\\tj \t J (k-l)\\ctj \t 



(-D'-'mV-.fe), (54) 



(k - 1)! \ct J \ct 



where U(X; c) is a rescaled version of U: 



u ^ c ^ cU ^ = xmr (55) 

Repeating the steps of this method with (|54[) . we get 

U k (t) = v fe • w fe! (56) 

where (vfc).; = (vfe)i/c l and is the solution to the matrix equation GfeW^ = e k , where {Gk)ji = (Gk)jiC 3 ~ l . 
To compute these products and ratios, a procedure as the one describes in the previous remark should be 
done to avoid overflow. We found the choice c = t^ 1 useful. 



2 As mentioned in I13| . the linear combination (1511 obtained by polynomial interpolation is equivalent to taking the linear 
combination which cancels the first n remainder terms in 1 148 I I. 
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3.3 Testing the methods 



In Tabic [TJ the Post-Widdcr method is tested in cases where U can be computed explicitly (see Section 0] 
for explanations of each example). The values in the table were obtained using a threshold of e = 10 -8 . 
With the exception of starred entry in Table [1] all converged within this threshold with n < 9. If sufficient 
convergence did not occur, Pg(t) is used as the approximation. In these cases, the approximation is very 
good. 

For the Poisson process, U is discontinuous at the integers. Therefore, we used instead t = 1.1, 10.1, 100.1. 
Since the Post-Widdcr method does not apply when U is discontinuous, we do not expect to get good results 
in the Poisson case. And indeed, the absolute errors are large. For instance, Table [2] shows that when 
t = 10.1, the algorithm converged but to the wrong value (the absolute error between the true and computed 
value is 0.4). The plot in Figure Q] shows in more detail the erratic behavior of the Post-Widder method in 
this case. Notice that for larger times this method approximates U (which in this case is a step function) 
with a straight line. 



Subordinator 


t = 0.01 


t = 0.1 


t = 1 


t = 10 


t = 100 


a-stable (a = 0.5) 


2.5 x 10" 13 


8.2 x 10~ 13 


1.7 x 10" 13 


7.4 x 10" 13 


2.6 x 10~ n 


Uniform Mixture of a-stable 


7.8 x 10" 12 


6.6 x 10~ n 


1.4 x 10" 10 


1.5 x 10" 9 


8.7 x 10" 10 


Gamma Process (k = 7 = 1) 


1.7 x 10~ 13 


7.0 x 10~ 14 


1.2 x 10~ 12 


3.4 x 10~ 10 


9.6 x 10~ 12 


Inverse Gaussian (7 = 6 = 1) 


1.3 x 10~ 5 * 


1.7 x 10" 12 


1.3 x 10~ n 


6.7 x 10- 11 


3.0 x 10- 9 



Table 1: Absolute errors between exact values of U and those given by this method for a selection of 
subordinators whose mean first passage time can be computed exactly. Here we used a threshold of e = 10 -8 . 



In the cases where U has little regularity, numerical integration (method 1) is the more accurate method. 
To see this, we applied the numerical integration method (using the MATLAB function quad) to the case 
of the Poisson process and obtained the absolute errors in Table [5] To compute these, we approximated 
the oscillatory integral in (|19p by partitioning the range of integration (0, 00) into intervals of the form 
Jo = (0, 7r/(2f), Ik = [kir/(2t), (k + 2)ir/(2t)}, k > 1, k odd, and then summing the integrals over each Ik 
until the contribution over one such interval became less then a threshold e > 0. Due to the slowly decaying 
nature of the integrand, obtaining convergence for e < 10~ 6 is difficult. 

We've found that integration is not feasible in many other cases (except possibly the a-stable case). 



Subordinator 


t = 0.01 


t = 0.1 


t = 1.1 


t = 10.1 


t = 100.1 


Poisson Process (integration) 


1.08 x 10~ 6 


1.57 x 10~ 6 


3.6 x 10~ n 


4.25 x 10~ 6 


1.09 x 10~ 4 


Poisson Process (Post-Widder) 


c i0" 16 


7.1 x 10" 10 


3.7 x 10" 2 * 


0.4 


0.4 



Table 2: Absolute errors between exact values of U and those given by numerical integration for the poisson 
process (whose mean first passage time can be computed exactly). Here we used a threshold of e = 10~ 6 . 
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Figure 1: A comparison of the true value of U for a Poisson process with no drift and its approximation 
given by the Post-Widder method. Since U in this case has little regularity, the quality of the approximation 
is not always good. 

4 Obtaining U(t) = KE(t) and Corr(£'(s), E(t)) for various inverse 
Levy subordinators 

We now present various examples of Levy subordinators {D(s), s > 0} and their inverses {E(t),t > 0} and 
calculate the one time moment U(t) = KE(t) and the correlation function cott(E(s), E(t)). We found the 
table of Laplace transforms [55] useful for some of the following calculations. 

Since we will encounter Laplace transforms which cannot be inverted analytically, we will study their 
asymptotics using a Tauberian Theorem ([9] page 10), which we state here for convenience. Recall that a 
function £(t), t > 0, is slowly varying at (respectively oo) if for all c > 0, lim(£ (ct) / £(t)) = 1 as t — > 
(respectively t — > oo). 

Theorem 4.1 (Tauberian Theorem) Let £ : (0, oo) — + (0, oo) be a slowly varying function at (respectively 
oo ) and let p > 0. Then for a function U : (0, oo) — > (0, oo), the following are equivalent: 

(i) U(x) ~ x p £(x)/T(l + p), x — > (respectively x — > ooj. 

(ii) U(X) - X-P-^il/X), A -> oo (respectively A -> 0J. 

4.1 Poisson process 

Consider the process -D(s) = /us + N(s), where {N(s), s > 0} is a Poisson process with rate r > 0. The Levy 
exponent for this is given by 

0(A) = p\ + r(l-e- x ), (57) 

and the Levy measure Tl(dx) is a point-mass at x = 1 with weight r. The sample paths of D are straight lines 
with slope /x, together with jumps of size 1 which happen at random times {rk}^ 1 , such that {t^ — Ti_i} < *L 1 
are iid exponential with mean 1/r. 

First consider the case with no drift, p = 0. Figure [2] displays a sample path of the Poisson process N(s) 
together with its inverse E(t) which is the first time N{s) exceeds the level t. Notice that each segment in 
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the plot of E has length 1 . For any < <o < 1 , one must wait a random time n for the process N to surpass 
to, thus E(to) = T\. More generally, the inverse subordinator E{i) is a sum of [t + lj iid exponential random 
variables with mean 1/r, implying E(t) ~ T(\t + lj , 1/r). 



Figure 2: A sample path of the Poisson process iV(s) together with its inverse. 
Using the density of a Gamma random variable, the 7-moment of the E(t), with 7 > 0, is given by 

Jt+lJ roc 

EE (ty = ^ru I 1 n / x^x^+^^e-^dx (58) 

(59) 
(60) 



r(L* + ij) Jo 

rL t + 1 Jr(Lt + lJ +7) 
r(Lt + lJ)rL*+iJ+7 

r(L* + i]+7) 

rTflt+lJ) ■ 



Setting 7 = 1 yields U(t) = [t + lj jr and the renewal measure, dU(t) , for the Poisson process is the measure 
which assigns a mass of 1/r to each integer n > 0. The covariance, (|17p . of this process is then given by 

CovOE( S ), = ^f: J + j - U(s)U(t). (61) 

Now assume positive drift, i.e. /i > 0. From (fTTj) and (|57|) . the Laplace transform of £/(£) is given by 

0W " At/tA + r(\ -fi~ A )) " < 62 > 

Due to the simple form of E/(t), can be calculated by numerical integration (i.e. method 1 in Section 
[3|) H. In Figure [3] we plot U(t) = EE(t) for various values of /1, namely fi = 0,0.1,0.5,1. We also compute 
the correlation coefficient. It is obtained using (|61[) in the case of no drift but it is not known in closed for 
fj, > 0. To compute for fj, > 0, we use (fTT)) . The integral in ([IT)) involves the renewal measure dU, but since 
in this case the function U(t) is strictly increasing and continuous for [i positive, dU is given by U'(t)dt. We 
obtain U'(t) here by discrete approximation. The correlation coefficient coir (E (s), E(t)) is plotted in Figure 
[3] as a function of t with s = 10 fixed for zero and non-zero drift. 

3 The Post-Widder method also works well when jj,/r 2> 0. It fails to converge however, when [i/r is close to 0. 
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Figure 3: (Left) Plots of KE(t), where E is the inverse of a Poisson process with rate r = 1 and various 
drifts. For fi > 0, these plot were generated using the numerical integration method in Section [3. II (Right) 
The correlation coefficient corr(I?(s), E(t)) of the Poisson process with s = 10 fixed and t varying with drifts 
fi = and fj, = 0.1. This was calculated using (|61| in the case \i = and by numerical approximation of the 
integral in equation |(T7|) for /j, > 0. Since here s = 10, the correlation equals 1 when t = 10 as well. 

4.2 Compound Poisson processes 

Let £i i = 1 . . . be positive iid random variables with probability measure v. The compound Poisson process, 
X(s), s > is defined by 

JV(s) 

D(s) = 6, (63) 

k=l 

where N(s) is a Poisson process with rate 1. In this case, the Levy measure II in is given by the 
distribution of the random variable £i. Therefore, the Levy exponent for this process is 

/•oc 

4>{X) = / (1- e- Xx )v{dx). (64) 

Notice that X(t) = for all t < Ti, where n has an exponential distribution with mean 1. This implies 
that £7(0) will also have an exponential distribution, and thus E(0) ^ a.s. The fact that E(0) ^ is 
characteristic of compound Poisson processes. Indeed, suppose that the inverse of a subordinator D satisfies 
E(0) > a.s., then in particular, EE(t) ^c>0at£^0 since E has cadlag paths. Thus, by the Tauberian 
Theorem, 

U{\) = ~ C[c]{\) = j, A^oo. (65) 

Thus, 4>(X)/c — > 1 as A — > oo, implying <fi is bounded, which means D is a Compound Poisson process (see 
[9], Corollary 1.1.3). 
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Now focus on the special case when £1 has a Pareto distribution, meaning £1 has probability density 



v(x) 



ax 



for x > 1, and a > is fixed. The Levy exponent is given by 



(i 



- \x 



)x a l dx = 1 — a Ei 1+a (A), 



(66) 



where Ei Q+ i(A) = ^ Q+ i dx is the Exponential integral. Using a series expansion^ of Ei Q +i, we have 



-aT(-a)\ a < a < 1 
0(A) ~ { (1 - 7e - log(A))A a = 1 
^rA a > 1 



A -> 0, 



(67) 



where j e w 0.5772 is the Euler constant. Since U(X) = (\<fi(\)) 1 , the Tauberian Theorem gives the large 
time behavior of the first passage time for this Compound Poisson process 

1 



U(t) 



-aT(-a)T(l + a) 
1 



t a , < a < 1 



1-7, 

a-1 
a 



>. + log(t) 
t, a > 1 



t a = 1 



(68) 



To calculate U, we found the Post-Widder method in Section[3]to be most usefuQ. In Figure|4]we plot U 
for a = 1/2, 1, 2, and in Figure [5] we plot U along with the asymptotic expressions given in (f68|) with log- log 
scale. The correlation corr(£'(t),^(s)) is plotted in Figurc[Hl Note that the renewal measure dU will give a 
mass of weight 1 to t = since U has a jump of size 1 there. For t > 0, (if/ is a measure which has density 
U'(t) which can be calculated as U is with the Post-Widder method as noted in a remark in Section [ 



4.3 "Mixture" of a-stable subordinators 

We consider here a continuous mixture of a-stable subordinators with < a < 1. Namely, the subordinator 
whose Levy exponent is given by 

0(A) = l\{P)\ p d(3 (69) 

- A:r 



(1 - e- xx )g p {x)dx, (70) 
where p is a probability density on (0, 1) and the density g p of the Levy measure is given by 

9v( x ) = J Y^W) md(3 (71) 

In general, this subordinator has no finite moments. The a-stable subordinator corresponds to the choice 
p(P', a) = 5{(3 — a) in (|69| . Here we will consider two extensions of this, namely when we choose a sum of 
two a-stables, p(/3; oc\, 02) = C\8{f3 — a.\) + C2<5(/3 — a%), with a-i < a\ and C\ + C2 = 1, as well as what we 
call the "uniform mix", which corresponds to the choice p((3) = 1 on (0, 1). 

4 obtained using Mathematica 

5 The only problem with this method occurs near t = 1 where this technique is slow to converge. The reason for this seems 
to be that U is not differentiable here. 
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Figure 4: The mean first passage time of the compound Poisson process with Pareto jumps computed using 
the Post-Widder inversion method for a = 0.5, 1, 2. 

4.3.1 Single a-stable 

Many properties of the inverse of an a-stable subordinator with < a < 1 are known (see for instance, [7]), 
but we restate them here for convenience. Consider the a-stable subordinator {D a (s), s > 0} with Levy 
exponent <f>(X) = X a . Taking the inverse Laplace transform of (fTTj) . we have that the mean first passage time 
for an inverse a-stable subordinator is given by 

1 (1 + a) 

which implies that the density of the renewal measure is U'(t) = t a ~ l /T(a). With this, the covariance can 
be given in closed form. Assume s < t, then 



^W.^W) - r( i + l)r ( a) £ '" - T »° + " - Tf ) ^ - WTa? < 73 > 
= r^ + s °'°(r(I^ F ( '- a ' a + l!s/< )-r(TT^)- < 74 > 

The above integral was computed in Mathematica and F denotes the regularized confluent hypergeometric 
function. 

4.3.2 Sum of two a-stable subordinators 

Here we consider the case when p{0) = C\S(P — ax) + C2^(/3 — 0.2) and the Levy exponent is given by 

4>{\) =CiA ai +C 2 \ a \ (75) 

where a.2 < ai and C\ +C2 = 1. This corresponds to the subordinator given by C\' /ai D ai (s) + C\' ,a2 D ' a2 (s) . 
From (fTT|) , the Laplace transform of the mean first passage time is given by 

^ (A) = ^A^^TT - (76) 
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10" 



-a = 0.5 
a= 1.0 
-a = 2.0 



t 



Figure 5: A log-log plot of the mean first passage time of the compound Poisson process with Pareto 
jumps computed using the Post-Widdcr inversion method for a = 0.5,1,2. The dotted lines represent the 
asymptotic expressions given in ([68]) . We see good agreement for large times. 



Using the Taubcrian Theorem, U has the following behavior in the limits t — > and t — > oo, 



r 



u(t) ~ { 



Cir(l + ai)' 
[ C 2 T(1 + a 2 ) ' 



t -> 



t — > oo. 



(77) 



Thus, we see a cross-over in power-law behavior (which is displayed with a log-log plot in Figure [7]). A closed 
form for the inverse Laplace transform of (|76j) could not be found, however U can be calculated numerically 
using the Post-Widder method (method 2 in Section [3]). Figure [7] shows plots of this function for various 
parameter values. 

To obtain the correlation, we set dU(t) = U'(t)dt, calculated U'(t) using the Post-Widdcr method and 
then evaluated the integrand in (|17p using numerical integration. The correlation is plotted in Figure [51 



4.3.3 Uniform mixture 

Here we consider the case when p(f3) = 1, for (3 G (0, 1). The Levy exponent is given by 

0(A)- f \Hfi= ±^r. 

Jo l °sW 



(78) 



The mean-first passage time for the inverse of this process has a closed form expression given by (see [25] 
Eq. 4.1.9) 



U(t) = cr 



c 



-i log(A) 
A 2 - A 



7e + e*r(0,i) + log(i). 



Here r(0,t) is the incomplete gamma function given by T(0,t) — Jl°° e z z 1 dz. Since, for t large, 

e'^-^z^dz < J e- z dz = 1, 



(79) 



(80) 
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Figure 6: The correlation coefficient con (E (t), E(s)) with s = 10 fixed for the compound process with 
Pareto jumps for a = 0.5, 1, 2. 



we obtain the "ultraslow" growth: 

U(t) ~ log(i), t^oo. 
The density of the renewal measure also has a nice form: 

U'{t) = e*r(0, t) - e^-H' 1 + 1- 1 = e*r(0, t). 



(81) 



(82) 



Using this, we can numerically calculate the covariance for the inverse of the uniform mixed subordinator 
with (jTTJ). The correlation is plotted in Figure [5] 

4.4 Generalized inverse Gaussian Levy processes 

The generalized inverse Gaussian (GIG) distribution [6], |llj is a distribution characterized by three param- 
eters (5, 7, k and has probability density function given by 



Here, K K denotes the modified Bessel function of the third kind ([H], section 3.2). The parameters of this 
distribution may take the following values: 




x > 0. 



(83) 



6 > 0, 7 > 0, if k > 0, 
5 > 0, 7 > 0, if k = 0, 
S > 0, 7 > 0, if k < 0. 



(84) 
(85) 
(86) 



Important subclasses in this family are 



• k>0, (5 = 0, 7>0 gives a Gamma distribution T(k, 2/~f 2 ) with density 



PGIG(0,~f, K )( X ) = 9n T ( K \ xK le -^ X i X > 0,K > 0,7 > 0. 




(87) 
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Figure 7: Plot of U(t) calculated with the Post-Widdcr method for a single a-stable, a sum of two a stables, 
and a uniform mixture. On the right we used a log-log plot and dotted lines to show asymptotic behavior 
given by (|77|) for the sum of two a-stable case. 

• k < 0, 5 > 0, 7 = gives a reciprocal Gamma distribution RT(k;, 6 2 /2) with density 

£-2re s2 

Pgig(sm,k)( x ) = 2-T(-K) xK ~ le ~^' x>0,k<0,5>0. (88) 
This distribution only has finite moments of order less than 

• k = — i, S > 0,7 > gives an Inverse Gaussian distribution IG(5, 7) with density 

p GIG{Stl ,- m {x)= 6 -^=x- 3 ' 2 e-^ 52x ~ 1+ ^ x \ x>0, 7 >0 ;( 5>0. (89) 
v Zir 



The Inverse Gaussian distribution is the distribution of the first time a Brownian motion with variance 
5 and drift 7 reaches the level 1 ([3], Example 1.3.21). All moments are finite if 7 > 0, but if 7 = 0, it 
becomes a totally right-skewed ^-stable distribution which only has finite moments of order less than 
1/2. 

The GIG distribution was shown to be infinitely divisible in [5] and its Levy-Khintchine representation 
was derived in Section 5 of [11] . The Levy-Khintchine representation given in [TT] is not in the form ^ and 
some computation is needed to bring it into this form. This is done in appendix IA1 

The Laplace transform of PgiG(S,j,k) f° r <5 > and 7 > is given by 

PGiGis^W = {^TT^x) KJh) (90) 

= exp(-</> GJG? ( 5j7iK )(A)), (91) 
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Figure 8: A plot of the correlation coefficient with s fixed at 10 and t varying for the a-stable process with 
a = 1/2, the sum of two a-stables with a± = .75 and 012 = -25 and the uniform mix. The correlation for the 
a-stable case can be calculated exactly using (|T4|) . In the other two cases, we numerically approximated the 
integral in (fT?"]) . 



where the Levy exponent 4>gig is given by 

/>oo 

;(<5,7,«)( A ) = / i 1 ~ e ~ Xx )9GiG(s.~f, K )(x)dx, (92) 



<J GIG( 



Thus, the Levy exponent is of the form ([2|) with drift fj, = 0. The corresponding Levy measure (density) is 
(see pj] or Appendix |A|) 



goicx^x) = — \J q ^jr^vm + Y ^5^)) dy + max(0 ' K) ) ' x>0 - (93) 

Here, J„ and Y v denote the Bessel function of the first and second kind, respectively, with index v ([14j. 
Chapter 2). 

By letting 8 — > with « > and 7 — > with k < in (|90|) . we obtain the Laplace transform of the 
gamma distribution (|87[) and the reciprocal gamma distribution (|88[) respectively. Doing so gives (see [11] ) 



PGJG(o, 7 ,«)(A) = ( 1+ yJ k, 7 >0, (94) 

i/2 



(?) ^ , 



PG/G(«,o,i.)(A) = r K- H (6V2\), k<0,5>0. (95) 

T(-K) 

The corresponding Levy measures are obtained similarly using ([93)) . (recall that Y v (z) — > —00 as z — > 0): 

ft t 2 

5g/g(o, 7 , k )0z) = - e ~~ x ( 96 ) 
1 f 00 e~ X3/ 
«./o ^ y(^| K |(oV2y) +^j«|(oV22/)) 
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The mean of a GIG distribution, when it exists, can be calculated from these Laplace transforms. If 
X ~ GIG(S, 7, k), then there are three cases for EX: 

EX = 5K X + (^8) kgRj7>m>0 (98) 
2k 

EX = — , (5 = 0,k>0,7>0 (99) 




-1< k < 0,5 > , 
EX = i i ' o " 100 

We talked so far about the GIG distribution. We now consider the corresponding Levy process, namely, 
the subordinator {-Dg/g(<5,7,k)( s )i s > 0} with Laplace transform 

Eexp(-A-D G / G ( 5i7)K )(s)) = exp(-s0 G/G (^ 7!(t) (A)). (101) 

From l|92p. the drift of this process is and its Levy measure is given by H(dx) = gGic(x)dx. Notice that 
this Levy process is indeed a subordinator since its Levy measure is concentrated on the positive axis and 
from (|M|) and (|9T|) . it follows that J Q (lAx)H(dx) < oo. Let Egig{5,^,k) be the inverse of D G iG(s,-y,K) and let 
UGiG{6,-y.K,)(t) = ^^GiG(5,-y, «)(*)■ The gamma process (5 = 0) and the inverse Gaussian process (k = —1/2) 
are treated in [21j . where closed form expressions are given for the renewal measure. In general, we cannot 
write ?7 G / G (t) in closed form, however, we can obtain expressions for large and small time behavior using 
the Tauberian theorem. 

Asymptotics of Ugig(s.i,k)i t "~ * 0: One has K K (x) ~ \f^ e ~ x as x ~~ * 00 ( EL section 5.4). We first 
consider the case when 5 > 0. When 7 > 0, (|9"0|) implies that as A — > 00, 




- log(K k (8^j 2 + 2\)) + log(K K (5j)) (103) 

= <V7 2 + 2A + 0(log(A)) (104) 
- 5V2A, A — > 00. (105) 

The case 7 = gives the same answer (by using (|95]l). If <5 = 0, we instead use (|94|) and get 

0GJG(O, 7 ,k)(A) - rclogfl + ^J (106) 
- Klog(A), A->oo. (107) 
Since Ugig(S,j,k)(^) = (A<^g/G(o,7,k)(A)) _1 , we get from the Tauberian theorem 

KlOg[t) 

t^0. (108) 

-4o y/t, s > 



UGIG(8,-y, K )(t) 



Asymptotics of Ugig{5,i,k)i t ^ 00: If 7 > 0, the renewal theorem implies that J7 G _r G (<5, 7 , re )(t) 
t/E£(l), where E-D(l) is finite and is given by 1(55] ) or ([99 ]) . 
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The case 7 = requires more work. We shall use the following series expansion of K K [x) as x — > with 
k > (see [23], page 121): 



f r(«) 

2 i-« 



aT" 1 



r(-«) 
4 K r(«) 



x jK +o(a; K ), 0<k<1 



1 



x _1 1 + -(2 7e - 1 + 21og(a?) - log(4))^ + o(x), k = 1 



T ^x- I 1 



(109) 



L2 1 -"' 



4(1 -k) 



X Z + o(aT K ), K > 1 



Now, (95]) and (fT09]) imply that for < -« < 1, 



<t>GIG(S,0, K )W = -log 

= -logf 

= -log 



2 4r 



.-2\ - K / 2 



r( 


-«) 


2 l+«/2 r 


-re A -re/2\ 


r(- 


«) J 


21+^/25- 


K A~ K / 2 \ 



if_«(5\/2A) 



log (if_ K (5\/2A) 



(110) 



(111) 



log 1 



V r(-«) 
r(«) 



log 



r(-«) 
2 i+ K 



5 « 2 «/2 A «/2 ) _ w ( x + 



r(«) 



2-«r(-«; 



2- K r(- K ; 



r(«) 



2-«5 2fi r(-K 

Similar calculations give for k = — 1, 

0G/G(5,O,k)M = 



A K , A^O. 



log (5V2A J pri((5V2A)j 
log (sV2X) - log(A'i(5\/2A)) 



(112) 

(113) 
(114) 



log (sV2\j - log (((5V2A)- 1 ) - log M + i(2 7e - 1 + 2 log(<h/2A) - log(4))(2<5 2 A) + o(x 2 ) 



--(2 7e - 1 + log(2<5 2 ) - log(4) + log(A))A, A - 0, 



(115) 



and for — k > 1, 



<I>GIG(S,0,k)W = -log 



-re/2 



A -re/2 



■log^ 
log 



r(-«) 

/ 2 1+k/2j-k^-k/2 

T(=k) 

/ 2 1+k/2 (5 -k A -k/2 



-if_re(5V2A) 



log A'_ K ((SV2A) 



(116) 



(117) 



V r(-«) 



log 



r ( K )c BOK / 



2 1+ K 



5 K 2 K / 2 A K / 2 ) - log ( 1 + 



4(1 + k) 



(2<5 2 A) + o(A) 



2(1 + AS 



■A, A^O. 



(118) 
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Thus, using the Tauberian theorem, we have the follow large time behavior of the mean first passage 
time of the GIG process 



t, n e K,7 > 0,6 > 0, 



2- k S 2k T(-k) _ k 



UGIG(S n ,K){t) ~ < 



r(«;)r(i - k) 



r K , 7 = 0, -1< k < 0,6 > 



t — ► oo 



(119) 



5 2 /2 



2 7e - 1 + log(2<5 2 ) - log(4) + log(t) 



t, 7 = 0,«=l,d>0 



2(-/t- 1) 
52 



i, 7 = 0, s < — 1, <5 > 



While (|108p and (|119|) give the asymptotic of f7, closed form expressions for Uqig are not known. Using 
our methodology, we can compute Uqig numerically. In Figurc[51 Ugig is plotted for three sets of parameter 
values. Figure [TO] shows a log-log plot which includes the asymptotic curves (|108[) and (| 1 1 0|) . Since this 
subordinator is strictly increasing, the renewal measure is given by dU(t) = U'(t)dt, and U'(t) can also be 
computed using the Post-Widdcr approach. The correlation corr (E (t), E(s)) is can also be calculated as it 
was in the other examples and is plotted in Figure [TT1 with s = 10 fixed. 



-8 = 2,y=3/4, k = 1 
8 = 1,7 = 0,K = -1/3 
■8 = 0, 7=3/4, k= 1/2 



52 

CD 

3 3- 




Figure 9: Plots of Ugig for various values of parameters S, 7 and k. Each of these plots were generated using 
the Post-Widdcr method with e = 10~ 5 . 
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Figure 10: Plots of Ugig on log-log scales. To dotted curves correspond to the asymptotic expressions given 
by (|108|) (left), namely t — > 0, and (| 1 19|) (right), namely, t — > oo. We see in these limits, the approximations 
are very good. 

5 Conclusion 

In this paper, we developed two numerical methods for calculating the function U(t) = ¥,E(t), where the 
process {E(t), t > 0} is the first hitting time of a Levy subordinator {D(s),s > 0} with Levy exponent 
<f> given by @. The function U has been shown to characterize all finite-dimensional distributions of the 
process E and is useful for calculating moments of E, for example, the covariance (see equation (|17p). 

The Laplace transform of U has a simple expression in terms of <fr, namely, U(X) = (X(j)(X))~ . Thus, 
calculating U involves computing the inverse Laplace transform of this function. The first method described 
in Section [3.11 computes the inverse of this Laplace transform by approximating the Bromwich integral given 
by (|18p . This integral can be computed by rewriting the integrand in terms of the real and imaginary parts 
of <j). We give explicit expressions for Re(</>) and Im(</>) in terms of the drift /i and the Levy measure II of 
the subordinator D. 

The second method described in Section [3 . 2 1 computes the inverse Laplace transform of U using the Post- 
Widder inversion formula given in (|25p . This formula is usually difficult to use because it requires evaluating 
derivatives of high orders. However, using our methods, can be calculated in a reasonable fashion. As 
with the integration case, all terms in this approximation are given in terms of only the drift /i and Levy 
measure LI corresponding to the subordinator D. We tested both of these methods in cases where U can be 
calculated exactly and obtained accurate approximations. 

As an application of our methods, we considered three families of Levy subordinators D and calculated 
the mean and correlation of their respective inverse subordinators E. The three families considered were 
(i) Poisson and Compound Poisson processes (ii) Continuous mixtures of a-stable subordinators and (iii) 
Generalized inverse Gaussian Levy processes. In each example, either the integration or Post-Widder method 
was useful for calculating the mean U(t) = EE(t). Once we computed U, the correlation function of E can 
be computed by numerically approximating the integral in (|17p . Along with these numerical approximations, 
we also gave in each case, asymptotic expressions for U. 
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Figure 11: Plots of the correlation coefficient p = con(E(t), E(s)) for the generalized inverse Gaussian 
process with various parameters values. Here, s = 5 is fixed and t varies 

6 Guide to Software 

We have developed a MATLAB software package which computes U for the examples above, as well as a 
program for a user defined example, which is available from the authors. The package includes 4 programs 
which calculate U(t): 

• The Poisson process(invert_poisson .m) 

• Compound Poisson process with pareto jumps (invert_pareto .m), 

• Sum of two a-stablc proccsscs(invert_sumas .m) 

• The generalized inverse Gaussian Levy process(invert_gig .m). 

The density of the renewal measure U'it) is calculated in all cases but the Poisson process. To use these 
functions, invoke MATLAB and add the inversesub directory to MATLAB's working path by typing 

> > addpath ( ' /yourpath/ inversesub ' ) 

Here, "yourpath" is the path in which the directory inversesub is located (for example, in Windows, this 
might look something like "C:/myhomedir/inversesub"). 

Table [3] defines the required inputs for each function including the method used (numerical integral or 
Post-Widder), the required parameters, and the outputs one obtains. In each case, the default tolerance is 
set to e = 10~ 6 . To change this, simply add an optional argument to each function which gives the desired 
tolerance. For example, typing 

>> invert_poisson(l . 2 , , 1) 

gives the mean first-hitting time U(t) = KE(t) for a Poisson process at time t = 1.2 with drift /i = 0a rate 
r = 1, and a tolerance of 10 -6 . Alternatively, one can type 

>> invert_poisson(1.2,0,l, .001) 

to compute U instead with a tolerance of 10~ 3 . If the requested tolerance cannot be met, the program will 
return a message saying so as well as a crude estimate of the error. 
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The programs using the Post-Widdcr method also returns an optional estimate for the derivative of U, 
U'(t). For example, typing 

>> [U DU] = invert _gig(l, 1,0,-1/2) 

assigns the value U(l) to U and U'(l) to DU where U corresponds to the inverse of the reciprocal gamma 
Levy process. Here t = 1 and 1, 0, —1 are parameter values (see Table[3]). 
Each program also accepts vector inputs for t. For instance, 

>> [U DU] = invert^gig([l 2 3], 1,0, -1/2) 

assigns the vector [U(l), U{2), U(3)} to the variable U and the vector [U'(l), U'(2), U'(3)] to the variable 
DU. 

The program "invert_empty.m" contains all the code of the previous examples with the piece which 
computes <fr- n > missing. To use this program, add a function call to line 30 of the code which computes the 
derivatives of cf> for your case. 



Function Name 


Method Used 


Number of Inputs 


Meaning & order of inputs 


Outputs 


invert _poisson 


Numerical Integration 


3 


t : time t in U(t) 
mu : Drift fi in (|57|) 
r : Rate r in l|57p 


U(t) 


invert _pareto 


Post-Widder 


2 


t : time t in U(t) 

a : Exponent a in (1061) 


U(t),U'(t) 


invert _sumas 


Post-Widdcr 


5 


t : time t in U{t) 
al : a\ in (|75|) 
a2 : a 2 in ^ 
cl : Ci in (J75]) 
c2 : C 2 in {75]) 


U(t),U'(t) 


invert _gig 


Post-Widdcr 


4 


t : time t in U{t) 
delta : S in (JH3D 
gamma : 7 in (1831) 
kappa : k in (1831) 


U(t),U'(t) 



Table 3: Information about the 4 functions included in the software package. The parameters in the fourth 
column should be entered in the order of top to bottom, for example, for the sum of two a-stable case, one 
would type invert_sumas (t , al , a2 , cl , c2) . 
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A Appendix: Levy-Khintchine form for GIG distributions 

Here, we show that the Levy exponent corresponding to a GIG distribution can be written in the form j2]) 
with drift /i = 0. This is done in the following proposition. 

Proposition A.l Let Pgig(8,i,k) denote the probability density of the GIG(S,^, k) distribution. Then the 
Laplace transform ofpciG is given by 

Pg/G(5, 7 ,«)(A) = exp(-^ G/G(5!7 ^ K) (A)), (120) 

where the Levy exponent 4>gig * s 

poo 

0G/G(d\ 7 ,K)( A ) = / i 1 ~ e ~ XX )9GIG{8 n , K ){x)dXi (121) 



with Levy measure (density) 

{e-^r* / f 00 e ~ x v \ 
/ „ . t9 . - — dy + max(0, k) , x > 0, 5 > 
x y ttW, 2 k |(^) + ^|(^)) ') . (122) 
-e-^ x x> 0,5 = 
x 

Proof. From [TT] Section 5.2, pgig is given in the form (|120p . but 4>gig is expressed in following alternative 
Levy representation which depends on the parameters S, 7, k, 

Ux 5K ^- K ^ + / (e^-l-»Ax)ff G /G(«, 7l «) i, 7 >0,Kel 
7^k(7<5) Jo 



-<t>GIG(6,~/,K,)W - < 



iX Tj2 + J ( e - 1 - iXx )9GiG(o,-fM x ) dx ' <$ = 0,k,7>0 

/•oo 1 _ -x 

i\S 2 / ff | K |(25 2 a;)dx 

iXx 



(123) 



+ / (e l x - 1 - iAa;l[ ,i](a;))s' G j G ( 5) o, re )(a:)da;, 7 = 0, n<0,8>0 



where, for v > 0, 



&>(z) = 2 rm ^ , V 2f / — \i > a; > 0. (124) 
tHzL/^U/x) + V(v^)] 

Thus, to show (|12ip . we much check that in each of the three cases in (|123|) . the drift term cancels with the 
"iAx" term in the integrand. 
Case 1: S, 7 > 0,k e M. 
We must show 

7^(7*) = io x 9GiG(S,f,n)dx. (125) 
Using (|122[) and performing the integration with respect to x gives 

*9Gi<K6„0** = max(0, K)j / o + ^ ^ - ^ ^ cfod* (126) 

2 f°° 1 

= — max(0, «) + / — T „ —^——^-dy. (127) 

7 2 lo ^( 2 / + 7 2 /2)y[J| 2 t |(v / ^)+>| 2 ; |(v^)] 
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We now use the change of variables y — > 26 2 y in the integral above and use the function g v defined in ()124[) 
to obtain 

f°° 1 f°° 1 

L ^ + 7V2MJ, K|2 (^) + ^ | (^)] dy ^V ^V 5 '^ (128) 
Now, we apply the integral representation given in formula (5.2) to obtain 



oo 



62 1 -^^ d y= SJ ¥wv- ^ 

Jo y + ^T lK\ K \[0V 
Thus, (|125p follows if we can now show 

— i z\ = T max (°> ) + — jj fx \ » K€lR,d,7>0. (130) 
7A K (7d) 7^ 7A| K |(d7) 

For this, we require the following two properties of the Bessel function K Ul (see for instance, [14] , formulas 
(3.15) and (3.22)) 

K v (x) = K- V (x), x>Q 7 veR (131) 
xK v+2 {x) = xK v {x) + 2(1 + v)K x + v {x) i>0,i/£l. (132) 

For k < 0, (|130p follows immediately from (|131[) . For k > 0, (| 1 32[) gives 

A maxM + fcg> = (13 3) 



7 2 ^ K (^7) 
7^ K (<57) 



This verifies (|125p . and hence finishes case 1. 

Case 2: S = 0,7, k > 0. 

This case is immediate, indeed, 



/ xg GIG{s ,o. K )(x)dx = 
Jo Jo Jo 



(134) 
(135) 



c\ roo / \ o \ poo 

i\— + / (e lXx - 1 - i\x) -e-^ x )dx= {e lXx - l)g GIG(0iJiK) (x)dx. (136) 
1 Jo \ x J Jo 



Case 3: 7 = 0, 5 > 0, n < 0. 

For this, we need to check that 

g\ K \(25 2 x)dx = / xg G iG(s,o,K)(x)dx (137) 

/o x Jo 

This follows by changing the order of integration and using the definition of g v : 

»1 roo p ~yx 

-dydx (138) 



Jo wy[J^{8^y)+Y^{S^m 

00 1 - e~ y 1 

y ipl+^M* (139) 

^ / — — 5|«|(2^ 2 y)^. (140) 

Jo y 
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This finishes the proof. 
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